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We unveil the role of substrate wettability on the reversal in the sign of the interfacial charge 
distribution in a nanochannel in presence of multivalent ions. In sharp contrast to the prevailing 
notion that hydrophobic interactions may trivially augment the effective surface charge, we demon¬ 
strate that the interplay between surface hydrophobicity and interfacial electrostatics may result in 
a decrease in the effective interfacial potential, and a consequent charge inversion over regimes of low 
surface charges. We also show that this phenomenon, in tandem with the interfacial hydrodynamics 
may non-trivially lead to either augmentation or attenuation or even reversal of the net streaming 
current, depending on the relevant physical scales involved. These results, supported by Molecular 
Dynamics simulations and experimental data, may bear far ranging consequences in understanding 
complex biophysical processes and designing nanofluidic devices and systems involving multivalent 
counterions. 

PACS numbers: 


A charged substrate in contact with a multivalent ionic 
solution may attract larger numbers of ions possessing 
opposite sign (counterions) than that necessary to screen 
its surface charge. This counterintuitive phenomenon, 
leading to a flip in the sign of the interfacial charge, is 
commonly known as charge inversion (Cl) [ij. Cl plays 
a vital role in biological processes such as DNA conden¬ 
sation 0, viral packaging Q, and drug delivery Q. Cl 
has also been found to be responsible for the reversal in 
the sign of the ionic current due to pressure-driven ad- 
vection (streaming currents) in the presence of an electri¬ 
cally charged interfacial layer adhering to the substrate 
Besides this. Cl is also responsible for an inverted 
mobility of charged species at high bulk ionic concentra¬ 
tions 0. 

Theoretically, the apparent anomalous behavior asso¬ 
ciated with Cl has been related to nontrivial ion-ion cor¬ 
relations, as verified through molecular scale simulations 
[M3- However, the reported theories on Cl fTll - [l3| have 
trivially considered the substrate to be perfectly wetting 
in nature, which is unlike the characteristic of real bio¬ 
logical and engineering systems. As a consequence, the 
role of substrate wettability on charge inversion remains 
poorly understood. Such inadequacies in theoretical per¬ 
ception possibly stem from the complexities in simulta¬ 
neously capturing the coupled electrostatics of Cl and 
the underlying wettability-induced interfacial phenom¬ 
ena over the relevant spatio-temporal scales. 

In this work, we bring out the inconspicuous role of 
wettability-modulated Cl on electrokinetic phenomena in 
nanofluidic channels. In contrast to the earlier reported 
findings on monovalent ions that hydrophobic interac¬ 
tions lead to an enhanced effective zeta potential [SO- 
we show that for multivalent ions, reduction in the sub¬ 
strate wettability may lead to a decreased effective zeta 
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potential, and even Cl at a relatively low surface charge 
density. In addition, we demonstrate that the interplay 
between hydrodynamics and interfacial electromechanics 
for multivalent counterions may non-trivially lead to ei¬ 
ther reduction (including inversion) or augmentation of 
the streaming current, depending on the physical scales 
involved. Our theoretical model devised on phase field 
formalism is shown to be in excellent agreement with 
molecular dynamics (MD) simulations as well with re¬ 
ported experimental data [^ , in an effort to support these 
unexpected trends. 

Towards establishing the interplay of hydrophobic in¬ 
teractions and Cl, we first refer to the schematic depicted 
in Fig. [TJ which qualitatively represents the solvent den¬ 
sity, ionic distribution and ionic transport characteristics 
in the near-wall region of a channel containing multi¬ 
valent ions. In our modeling paradigm, we couple the 
electromechanics with the physics of hydrophobic inter¬ 
actions over interfacial scales by relating the potential 
distribution (ip) in the wall-adjacent layer with the sol¬ 
vent density distribution through the Poisson equation: 
^ = -J2zien°goi{y), fory > a/2, where y 

is the wall-normal coordinate, a is the diameter of hy¬ 
drated ions (assumed same for both the ions), e is pro¬ 
tonic charge, Zi and are the valency and bulk number 
density of ions of species, goi (y) is a modified wall- 
ion distribution function for the species (see [13 for a 
discussion on the various interactions affecting gcu); e is 
the permittivity of the medium which varies spatially in 
accordance to the relative phase distribution originating 
out of hydrophobic interactions [l^l- The near-wall sol¬ 
vent density variations due to hydrophobic interactions, 
steric effects and ion-ion correlations [13 are simulta¬ 
neously accounted by a modification of the Boltzmann 
description of the number density of the ions , so that: 

goiiy) = ^Ciiy) exp {-PzieLiif) - (F - Fo)) , 

where L^ip) and F are functions of potential (ip) and 
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FIG. 1. (color online) Schematic depicting the two scenarios: 
the wall is (a) non-wetting and (b) perfectly wetting. For 
case (a), the solvent depletion (represented by the gradient in 
the background color) leads to a local variation in the solvent 
density (leftmost plot; the dashed line depicts the non-wetting 
case where there is a depleted wall-adhering phase due to hy¬ 
drophobic interactions). The consequent change in the ionic 
distribution, (gop) is shown in the rightmost plot. Owing to 
lower solvent density, there is a larger counterion concentra¬ 
tion closer to the wall for case (a), denoted by the dashed 
line. Cl is more strongly manifested in case (a) as compared 
to case (b). Notably, over length scales of few nanometers, a 
larger ionic flux in presence of an external pressure-gradient 
driven flow is also obtained for case (a), owing to a less-dense 
phase which leads to slipping hydrodynamics over interfacial 
scales with excess mobile charges. 


transverse coordinate respectively, Q {y) is the exclusion 
volume term of the ion, whereas Fq = 1/(1 -|- kqu), 
and for a binary electrolyte solution, kq = lim k = 

y—^oc 

\/ ” ^efkbT^"^^ ' Here k is the local Debye-Huckel param¬ 
eter given by, = [vPe^/ {e{y) hT)) {Zpgop - Zngon), 
P = {kbT)~^, kb being the Boltzmann constant, T be¬ 
ing the absolute temperature, and = ZpTip = — 
where rip, n° and Zp, Zn are the bulk number densities 
and valency of cations and anions respectively [ni,p(2/) 
is the near wall solvent density distribution, and pi is 
the bulk solvent density. The Poisson equation, as men¬ 
tioned above, is solved in conjunction with the equation 
^ (e {y) = 0 for y < a/2 with the following match¬ 

ing boundary conditions for obtaining ipiy)'- (i) continu¬ 
ity of ijj and at y = a/2, (ii) known surface charge 
density at the wall, a = —{e (y) dV'/rf?/)ly=o sym¬ 

metry boundary condition at the centerline of the chan¬ 
nel, {dtp/dy)\y^f,^ — 0, where h is the half channel height. 

In order to mathematically close the above set of 
equations, we relate the variations in e (y) and p (y) /pi 
with the distribution of an order parameter (p, by ap¬ 
pealing to the phase field formalism [I^. Here (p = 
(n„ —ni)/{ny + n/), where rii represents the number den¬ 
sity of the phase [l^. cp = —1 represents the bulk 
liquid phase (denoted by subscript 1), and cp = 1 (de¬ 
noted by subscript v) represents the wall-adhering low 
density phase formed due to hydrophobic effects [^ . 
The calculation of equilibrium (p starts by first consider¬ 
ing a free-energy functional which represents the excess 
Ginzburg-Landau-like free energy for a binary mixture 


given by [2l|: AH ((/) = f 


i(g)’ + *“W 


dy + Hs, 


where Alo{(P) is the bulk free energy having a double well 
potential profile [13, [13| , Aw = A ^^2 _ ^ and Hg is 

the surface energy that takes into account the interac¬ 
tions between the substrate and the fluid M- Here B is 
a positive constant such that B ^ kbTc with Tq being 
the critical temperature for the liquid-vapor coexistence, 
whereas ^{d(p/dy)^ is the interfacial energy with k being 
a positive constant [2ll| . The minimization of free en¬ 
ergy functional results in the Euler-Lagrange equation: 

= 0 [131. The interfacial 
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value of (p is related to the surface wettability through 
the contact angle, which is given by cos^u, = 

[13 m. Combining the double well potential for Aw 
with the Euler-Lagrange equation yields the following 
governing differential equation for the order parameter 
variable in a dimensionless form: — C^{(p^ — l)(p = 0, 

where C = aj^ {k ^ B^^; ^ being the multiple of inter¬ 
facial thickness 0), and y = y/a. The correspond¬ 
ing boundary conditions are: {d(p/dy)\y^f^^^ = 0 and 
'P\y=o = 4>S- Properties, such as p and e are interpo¬ 


lated as: X = Xv^-^ + where x is a generic prop¬ 

erty. The value of parameters and properties used in the 
numerical simulation are Zp : z„ = 2 : —1, £/ = 78.5£oj 
where £o, is the permittivity of the vacuum, £„/£/ = 1/15 
1^,^/pi = 2.306 X 10-^ T = 298 K, C = 0.38 nm 
20L l25l |. The values of the remaining parameters are 
stated later. 

When the ionic solution, energetically described as 
above, is driven by an external pressure gradient applied 
orthogonal to the direction of the variation of (p and tp, 
one may have an effective axial electrical body force even 
though no external axial electric field is applied. This 
body force is because of the development of a back elec¬ 
tromotive potential (also known as streaming potential) 
in an otherwise pressure-driven flow field, and may be de¬ 
scribed as: Fek = PeEs, where pe is the electrical charge 
density and Eg is the streaming electric field. This un¬ 
determined streaming field is obtained by setting the net 
ionic current, which is the sum of the streaming (advec- 
tion) current and electromigration current, to be zero in 
presence of an applied pressure gradient [^, [^ . 

In Fig. 2(a), we depict the variation in the effective sur¬ 
face potential (tp = Petp') | __g ^ vs dimensionless surface 


charge density, cr = aaP je, for a 2 : 1 electrolyte solution 
at various concentrations and contact angles (see [13 for 
a detailed description on the theoretical treatment of the 
modified Poisson-Boltzmann formalism). Additionally, 
in Fig. 2(b), we plot wall-ion distribution functions for 
2 : 1 electrolyte solutions at contact angles of 0°, 40° and 
100°. We compare these results against reported Monte 
Carlo (MC) simulations for a perfectly wetting substrate 
0. Considering reported experiments on nanochannel 
made of fused silica Q, we have chosen the contact an¬ 
gle, depending on the relative fraction of Si-OH and Si- 
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FIG. 2. (color online), (a) Reduced potential at location 
‘a/2’, ^(1/2) = l3eil}{a/2) as a function of reduced surface 
charge density, a = aa?je for a 2:1 binary electrolyte at 
various electrolyte concentrations and contact angles, (b) 
Singlet wall-ion density distribution functions for 2:1 binary 
electrolyte with a concentration of IM and a — —0.1685 at 
different contact angles. Comparison is made against MC re¬ 
sults reported in (open circles). Inset depicts the counte¬ 
rion concentration spatial variation obtained from the present 
model and MD simulations (see for details pertaining to 
MD simulations). It also depicts a concentration threshold 
(same as maximum counterion concentration in MD), above 
which the ions are considered to be jammed and immobile. 
Here, Co-I and Coun-I represent the coion and counterion 
wall-ion distribution functions respectively. The hydrated 
ionic diameter is considered as 0.3 nm according to the MC 
simulations [3]. 


0-Si groups on the surface dictated by the fabrication 
technique, to vary between 0° and 40° [111. In addition 
to these contact angles, we also show the plots obtained 
for 9w = 100°, considering the emerging trends of using 
highly hydrophobic polymeric substrates for fabricating 
miniaturized fluidic channels [^ . 

We refer our results in perspective of the potential at 
the plane located at ‘a/2’(shear plane), also known as 
the zeta potential [2^ . For a perfectly wetting substrate 
{9w = 0°), the zeta potential calculated using the th eory 
is in good agreement with MC simulation results [ 3 , . 

When the counterions are monovalent, (positive half of 
Fig. 2(a)), two trends of '!/j(l/2) vs a are evident for an 
increasing contact angle. First, for low (t(0 < a < 0.1), 


the zeta potential increases with the contact angle; a phe¬ 
nomenon previously noted as well [HO- Second, for 
high ct(0.1 < a < 0.3), the zeta potential decreases with 
the contact angle as attributable to Cl at the high sur¬ 
face charge densities. However, when the counterions are 
divalent (negative half of Fig. 2(a)), the magnitude of 
zeta potential decreases with the contact angle, eventu¬ 
ally flipping its sign, after which the magnitude of this 
inverted zeta potential increases with the contact angle. 
All these observations indicate CL More importantly, the 
substrate wettability, in conjunction with the ionic va¬ 
lency and substrate charge density, plays an important 
role in modulating the CL 

For low surface charge densities, low ionic concentra¬ 
tions, and monovalent counterions, Boltzmann distribu¬ 
tion of ions is valid, implying no CL However, in the 
other spectrum of high surface charge density or multi¬ 
valent counterions, due to increased electrostatic inter¬ 
actions between substrate and counterions, a high coun¬ 
terion density exists near the substrate which can lead 
to CL Lower substrate wettability also increases counte¬ 
rion density near the substrate, thereby leading to Cl due 
to several reasons - (i) counterions can occupy the vacant 
sites left by the solvent molecules, (ii) reduced possibility 
of existence of hydrated ion near the substrate increases 
the attractive force between substrate and the counteri¬ 
ons, and (iii) high electric field near the substrate due 
to the reduced solvent permittivity further increases the 
corresponding attractive force. In conclusion, an inverted 
interfacial potential, which is due to Cl, can be observed 
even at low surface charge density provided the coun¬ 
terions are multivalent and substrate is hydrophobic as 
shown in Fig. 2(a) for the case of (0^ = 100°). 

Table U depicts the streaming currents obtained from 
MD simulation (see [13 for the details of MD simula¬ 
tions and streaming current calculations) of electroki- 
netic transport of CaCl 2 and MgCl 2 solutions in a chan¬ 
nel of 8 nm height, as compared against those obtained 
using the present model, for different contact angles. In 
this context, it may be noted that Cl has been experimen¬ 
tally observed for MgCl 2 and CaCl 2 solutions by stream¬ 
ing current measurements Q where a negative advective 
current beyond a threshold ionic concentration indicates 
Cl, albeit in larger channels (490 nm) that may be ex¬ 
perimentally probed unlike the ones investigated through 
MD simulations (8 nm). Continuum-based explanations 
of these experimental observations, however, have over¬ 
looked the effect of the substrate wettability on the per¬ 
tinent observations [1) O) Hi- Towards this, we obtain 
the streaming currents at the contact angles of 0 ° and 
40°, in accordance with the contact angles on fused sil¬ 
ica substrates [l^l, which are shown in Fig. [3] along with 
reported experimental data (denoted by open circles) Q . 

It can be seen from Table U that the streaming currents 
obtained using the present model are in good agreement 
with the MD simulation results. These results upscale 
nicely for relatively larger channels (490 nm) as well, 
and match well with experimental trends 0 (see Fig. [3]). 
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TABLE I. Variation of streaming currents per unit width (from both MD and present model) with the concentration of CaCL 
and MgCL solutions at contact angle of 40° and 110° for a nanochannel of height 8 nm and a — —0.15 C/m^.The hydrated 
diameters of Mg^"*" and Ca^"*" ions are considered as 0.6 nm and 0.52 nm respectively in accordance with the literature |30l. 1^. 


CaCb, a = 0.52 nm 



MgCb, 

a = 0.6 nm 



= 40° 

dw 

= 110° 


= 40° 

= 110° 

C 

Imd 

Icont 

% 

Imd 

^Cont 

% 

Imd 

^Cont 

% 

Imd 

^Cont 

% 

(M) 

(A/m) 

(A/m) 

error 

(A/m) 

(A/m) 

error 

(A/m) 

(A/m) 

error 

(A/m) 

(A/m) 

error 

0.5 

5.48 

5.46 

0.36 

88.58 

91.18 

2.94 

5.54 

5.4 

2.53 

88.58 

91.11 

2.86 

0.6 

5.42 

5.43 

0.18 

88.69 

91.16 

2.78 

5.43 

5.37 

1.1 

88.64 

91.08 

2.75 

0.7 

5.43 

5.4 

0.55 

88.74 

91.15 

2.72 

5.48 

5.34 

2.55 

88.73 

91.06 

2.62 

0.8 

5.4 

5.38 

0.37 

88.81 

91.13 

2.61 

5.45 

5.32 

2.38 

88.82 

91.04 

2.5 

0.9 

5.42 

5.37 

0.92 

88.83 

91.12 

2.58 

5.45 

5.3 

2.75 

88.85 

90.51 

1.87 

1 

5.44 

5.35 

1.65 

88.95 

91.1 

2.42 

5.43 

5.26 

3.13 

88.93 

90.47 

1.73 



FIG. 3. (color online) Variation of streaming current per unit 
pressure drop with the concentration of MgCL solution and 
for different contact angles in a channel of height 490 nm 
and a = —0.15C/m^. Solid, dashed and dash-dotted lines 
represent results obtained for contact angles of 0°, 40° and 
100° respectively. Symbols represent the experimental results 
as reported by van der Heyden et al. . Dotted lines represent 
the results obtained via the Poisson-Boltzmann route. The 
hydrated ionic diameter of Mg^^ is taken from |^. [^. 


As seen from the MD simulations and continuum pre¬ 
dictions in Table HI streaming currents increase with the 
contact angle for 2h = 8 nm. This is due to an increase 
in both the solvent velocity (see [13 for MD evidence 
of hydrophobicity—induced slipping hydrodynamics) and 
counterion charge density near the wall with a decrease 
in the substrate wettability. 

We note, however, that the streaming currents mea¬ 
sured in experiments and obtained using continuum pre¬ 
dictions reduce with the contact angle as shown in 
Fig- El This is because, at experimental length scales, 
the slip length is much smaller than the channel height 
(2h = 490nm), thus rendering the condition at the wall 
to be effectively no-slip [s^- Moreover, a layer of fluid 
near the wall could be immobile due to the increased vis¬ 
cosity of the solvent which itself is due to the jamming of 
the counterions [s^ . Therefore, despite an increase in the 


counterion charge density near the wall with an increase 
in the contact angle, a simultaneous increase in the width 
of immobile region near the wall due to crowding of coun¬ 
terions leads to a decrease in the streaming current with 
an increase in the contact angle over the reported experi¬ 
mental scales. In this work, the jamming layer is defined 
as distance over which ionic concentration is greater than 
a threshold concentration as depicted in the inset of Fig. 
2(b). Quite naturally, the results obtained from the PB 
equation show gross overestimations in streaming current 
as compared to the reported experimental data. Even 
though the streaming currents obtained by using a con¬ 
tact angle of 0° and 40° are close to the experimental 
results, the most important artefact of the experimen¬ 
tal observations — negative streaming currents at high 
concentrations, are never observed for a contact angle of 
0°. On the other hand, such negative streaming currents 
are observed by using a contact angle of 40°, thereby 
suggesting the focal role of substrate wettability towards 
consistent predictions of experimental observations. It 
can also be observed that the reversal of streaming cur¬ 
rents occurs at progressively lower ionic concentrations 
as the contact angle increases. Similar to the influence of 
substrate hydrophobicity, with an increase of ionic con¬ 
centration, more number of counterions get packed near 
the surface (as manifested through the decrease in De¬ 
bye length in the mean field theory) and this can lead to 
CL As the inversion of streaming current is a manifesta¬ 
tion of Cl, such an inverted streaming current can also 
be observed at lower ionic concentrations provided that 
the wettability of the substrate is less than a threshold. 
This is evident from Fig. |3l where the threshold ionic 
concentration for the inversion of streaming current at a 
contact of 100° (Ccr, 2 ) is lower than the threshold ionic 
concentration at a contact angle of 40° (Ccj-,i). 

To summarize, we have provided a new theoretical 
framework for revealing the role of substrate wettability 
on the electrokinetic transport through nanochannels, for 
a solution containing multivalent ions. We demonstrate 
the dual role played by substrate wettability towards af¬ 
fecting the electrokinetics — on one hand it leads to an 
alteration in the ionic charge density near the wall due 
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to solvent-substrate-ion interactions, while on the other 
hand it leads to a change in the current flux depending on 
the interplay of slip length and the characteristic length- 
scale of the channel. The streaming currents evaluated 
from such a framework are in excellent agreements with 
simulations and experiments, yielding a broader perspec¬ 
tive in the inconspicuous yet decisive role of the substrate 
towards the electrokinetic phenomena and biochemical 
ion transport processes prevalent in nature and engineer¬ 
ing. 

We would like to acknowledge Mr. Chirodeep Bakli for 
his help in carrying out MD simulations and insightful 
discussions regarding the same. 


Appendix A: Modifications in electrostatic potential 
beyond the Poisson-Boltzmann picture 

We first extend the traditional Boltzmann distribution 
of ions by accounting for the finite size of ions through 
an excluded volume contribution and the electrostatic 
correlations among the ions through a fluctuation poten¬ 
tial [Hill. The key assumption behind the modeling 
of these terms, as given in the literature [ll|i El; is to 
regard the solvent as continuous uniform dielectric with 
a bulk dielectric constant. However, this needs to be 
supplemented with additional artefacts because of hy¬ 
drophobic interactions. In particular, the loss of hydro¬ 
gen bonds near an extended hydrophobic surface causes 
water to move away from these surfaces, thereby pro¬ 
ducing a wall-adjacent depletion layer [^. According 
to Joly et al. El; the depletion of solvent near a hy¬ 
drophobic surface gives rise to an excess effective inter¬ 
action potential in addition to the electrostatic potential 
as Vext = —kbT\xi[p{y) /pi), where kj, is the Boltzmann 
constant, T(= 298 K in this case) is the absolute temper¬ 
ature of the solvent, p {y) is the near wall solvent density 
distribution, and pi is the bulk solvent density. This ex¬ 
cess potential contribution is calculated in the present 
theory by expressing p as a function of the order param¬ 
eter, and appealing to the numerically obtained profile 
variation of 4>. 


Appendix B: Governing equation for EDL potential 
distribution 


For a binary electrolyte {zp : z„), where the subscripts 
p and n denote the positive and negative ions respec¬ 
tively, the non-dimensional governing differential equa¬ 
tions and boundary conditions are obtained using y = 
yja, V) = f3etp, Li = f3eL, ipy = Bq = Bofie/a, 

4>f = and X = x/xi as the non-dimensional vari¬ 

ables and introducing the dimensionless parameters as 
yo = koa, r{y) = I3e^ /(4:TTe(p)a), where T (y) is the 
Coulomb coupling parameter [lll |. Due to finite ionic 
size, no charge exists within a distance of radius of the 
ion from the wall (Stern layer). Accordingly, the dimen¬ 
sionless Poisson equation reduces to m 


d 

dy 



(Bl) 


in the Stern layer, whereas in the diffuse layer, we have 

m 

= - — — [gop-gon], for y> 1/2 (B2) 

dy \ dy J Zp - Zn 

The corresponding wall-ion distribution functions are 
given as [ill 


gop = -C(y) 

Pi 


-ZpLiiPj) --—(F - Fo) 


(B3a) 


p{y) 

gon = - C{y) 

Pi 


-z„Li(^)-^^(F-Fo) 


(B3b) 


where El El; 


LM) = + 1) + ^ [1 + (^ - 1)(2/ - 1/2)] V’(O) 

-hi [(3/2 - y)(y - 1/2) + F[f - 5/4)] ^'(0) 

y-ti 

- ^^2 1/2 < y < 3/2(B4a) 

1/2 


y+i 

LiW = Y [i^iy +1) + i^iy - 1)] ~ 2 / '^^y'^^y' 

y-i 

for y > 3/2(B4b) 


F = {4, + /(53)/[4 -k Ka{l + 2y) + 

forl/2<y<3/2 (B5a) 

F = {1 + /d 2 )/[(l + Ka){l - /<5i)],for y > 3/^B5b) 

§2 = exp[2Ka(l — y)] sinh(Ka)/ (2Kay) (B6a) 


di = (52[Kacosh(Ka) — sinh(Ka)]/[(l -I- ko) sinh(«;a)] 

(B6b) 

(1-exp{Ka[(l-k 2y)^/^ - 1 - 2y]}) 

03 — -z- 

nay 

(l-2y+(l + 2y)'/2) 

-k-1-;-(B6c) 

y 


<54 = <53-(l/y)(l-kexp{Ko[(l -k 2y)^/^-l-2y]}) (B6d) 


0, for 0 < y < 1/2, 
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V 

C,{y) = H{y- i) exp i -2-Ka^ J dyi 

I h/a 


liY -yi)9op(X)exp[-Zp(l}f{yi,Y)]dY 
max(l/2,yi-l) 

yi+L 

+n° J{Y - yi)gon{Y)exp[-Zn(j)fiyi,Y)]dY 

\ max(l/2,yi-l) 


/ 


(B7) 


0 /(yi>^) = ■Bo(?/i) < exp(-/ta) + 


/ 


(1 + 4^1?) 


1/2 


exp[—Ka(l + ‘iyiY)^^'^] 


(B8) 


where f = (e — Sw) / (^ + Sw)-, £w being the permittivity of the wall. In this study, we assume, ew = e, so that 

/ = o pill. 


Bo (y) = {iPv/tt) 


2 [(1 + Kay + Ka) exp {—na) + exp (—Kay)] — 

{f /nay) [2Kayexp {—nay) + (1 + na) exp {—na) (exp {—2Kay) — 1)] 


,for 1/2 <y < ](B9a) 


Bo( 2 /) = {exp(Ka)/[47r(l +/ta)(l -/(5ii)]} for y > 1 


(B9b) 


(2y - imi/2) + {y + l/2)(y - 3/2)iA'(l/2) 

y+a _ 

+2'ip{y + l)-2 f 'il){y')dy' 

1/2 


for 1/2 <y <3/2 


(BlOa) 


Ipy = 2TT 


y+1 


i’iy +1) + ip{y - 1) - / '>P{y')dy' 


y-i 


for y > 3/2 


(BlOb) 


<5ii = 


(ko cosh (ku) — sinh (ku)) exp [ku (1 — 2y)] sinh (ko) 
2Kya (1 + Ka) sinh (kg) 


(Bll) 


The boundary conditions used to solve non-dimensional 
governing differential equations for i/ with y S [0,/i/a] 
are : (i) continuity of //> and at y = 1/2, (ii) 

known surface charge density at the wall, aea/ {siksT) = 
— (e (y) di/i/dy) | and (hi) symmetry boundary condi¬ 

tion at the centerline of the channel, {d'ip/dy) |= 0, 
where h is the half channel height. We reiterate here that 
the e appearing above is not constant; rather, it is de¬ 
termined by locally interpolating between the vapor and 
liquid permittivity as determined by the spatially vary¬ 
ing order parameter. Mathematically, e = 

These coupled, non-dimensional governing differential 
equations are numerically solved using commercial soft¬ 
ware COMSOL Multiphysics 4.4. 


Appendix C: MD Simulation Details 

The model system for MD simulation consists of a 
parallel plate geometry with a channel height of 8 nm 
and planar wall dimensions of 5 nmx5 nm. Periodic 


boundary conditions are applied in X, Y and Z; X be¬ 
ing the axial direction. The channel is filled with 6936 
water molecules (SPC/E model) and the number of ions 
is determined by the bulk concentration and overall elec¬ 
troneutrality of the system. After energy minimization, 
the system was equilibrated for 12 ns, followed by non¬ 
equilibrium simulation for another 12 ns. All the simu¬ 
lations are performed with time step of 1 fs. The wall 
molecules are kept fixed at their respective lattice posi¬ 
tions. The fluid (water and ions) is actuated by applying 
a uniform acceleration in the X direction. Wall wettabil¬ 
ity is altered by tuning the Lennard-Jones (LJ) parame¬ 
ters of the hetero-nuclear potential between the wall and 
the oxygen of water molecule. The temperature is kept 
constant at 300 K by a Noose-Hoover thermostat and 
the long range electrostatic interactions are calculated by 
Particle Mesh Ewald [s^. The ion distribution and ve¬ 
locity profiles are obtained by suitable binning and the 
streaming current is calculated using the methodology 
discussed in section VI. 
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Appendix D: Evidence of Slip from MD based 
velocity data 



FIG. 4. (color online) Mean and standard deviation of the 
molecular velocity profiles (symbols) in a channel of height 8 
nm and surface charge density of -0.15 C/m^, obtained from 
the MD simulation data, are shown at contact angles of 40°, 
80°, and 110°. The dashed lines represent the fitted Poiseuille 
flow profiles at the corresponding contact angles. Note that 
the deviation from the mean velocity at the contact angles 
of 40° and 80° is imperceptibly small of the fitted Poiseuille 
profile. 

Fig. m depicts the molecular velocity data obtained in 
a channel of height 8 nm, with surface charge density of 
-0.15 C/m^, at contact angles of 40°, 80° and 110° re¬ 
spectively. This MD data is matched with a Poiseuille 
flow velocity profile for which the pressure gradient and 
slip length is obtained from the MD data. The solvent 
viscosity is obtained by matching the velocity gradient at 
the wall for both MD and continuum predictions, as was 
also done in the literature [illii- Navier slip bound¬ 
ary condition at the channel walls, 

is utilized for obtaining the corresponding Poiseuille flow 
velocity profile, where Ig is the slip length. As seen from 
Fig.lH both the slip velocity and slip length increase with 
the contact angle. Since for a given wettability, slip phe¬ 
nomenon appears to be more pronounced for channels 
with reduced dimensions, one can safely remark that the 
maximum slip length for the investigated experimental 
scenario {2h ss 490 nm) must be lower than the same 
for the investigated MD scenario {2h Ki 8 nm). As the 
maximum slip length at the highest contact angle investi¬ 
gated is much smaller than the characteristic channel di¬ 
mension over the probed experimental scales (s» 490 nm), 
the no slip boundary condition remains justified over the 
reported experimental scales, despite the emergence of 
slipping hydrodynamics for the corresponding physical 
scales probed through MD simulations (s» 8 nm) 


Appendix E: Near-wall density fluctuations from 
MD data 



FIG. 5. (color online) Normalized number density of water 
molecules across the half channel height of a channel of height 
8 nm, surface charge density of -0.15 C/m^, at different con¬ 
tact angles. The circles, diamonds and square symbols repre¬ 
sent the data at contact angles of 40°, 80° and 110° respec¬ 
tively. 

In order to elucidate the influence of substrate wetta¬ 
bility on the slip velocity or the velocity profile, we plot 
the fluid density profile across the half channel height. 
Fig. [5] shows the plot of the normalized density of water 
molecules against the channel height at contact angles 
of 40°, 80° and 110°. For low values of contact angles 
{dw = 40°), the number density of water molecules ex¬ 
hibits a hump near the surface, which progressively gets 
attenuated due to the bulk thermal motion as seen from 
Fig. [5j At contact angles of 110°, the surface becomes 
more hydrophobic, so that the water density near the sur¬ 
face becomes smaller than that at lower contact angles 
as shown in Fig. [5] This decreasing near-wall water den¬ 
sity hallmarks the increasing interfacial slip as depicted 
in Fig. m 


Appendix F: Calculation of Streaming Currents 

Substrate wettability not only affects the electrostat¬ 
ics but also significantly influences the associated hydro¬ 
dynamics, explicitly through the density variations and 
implicitly through its effects on electrostatics, as exper¬ 
imentally observed through streaming current. Stream¬ 
ing current is defined as the advective pressure driven 
current of the mobile ions under the constraint of a net 
zero current (sum of advection current and electromigra¬ 
tion current) in the system [^ . The streaming current 

2hja 

is thus given hy I = Wa f pe (y) u (y) dy^ where W is 

0 
























the channel width, (y) is the local charge density and 
u (y) is the fluid velocity profile. Over length scales in the 
purview of MD simulations, the velocity profiles are cal¬ 
culated as explained in the Section IV. At experimental 
length scale of 2h = 490 nm, the slip length being much 
smaller than the channel height {2h = 490nm), the no 
slip boundary condition at the wall may still be consid¬ 
ered to be valid [S^. Moreover, a layer of fluid near the 
wall is likely to be immobile due to the increased vis¬ 
cosity of the solvent which itself is due to the jamming 
of the counterions^a phenomenon termed as Charge In¬ 
duced Thickening [33|. Therefore, over such experimental 


length scales, the velocity profile is given by, 

{ 0, for 0 < y < 7 

7? w - (S - 7) I) + J it‘‘I- 

for 7 < y < h/a 

(FI) 

where 7 is the non-dimensional distance over which the 
concentration exceeds the threshold concentration as de¬ 
fined in the main text (see the inset of Fig. 2(b) in the 
main text for the details about calculation of threshold 
concentration). Here the first term represents the flow 
due to applied pressure gradient, —dP/dx. The second 
term represents the flow due to induced streaming field 
{Eg) which satisfies an overall current electro neutrality 


Pet = aert 


n2hja n2hja 

/ (5op - gon) u (y) dy + {EgaX/2) / {gop + gon)dy = 0 

J J 0 


(F2) 


Streaming Current 


Conduction Current 


leading to 


a'^en^ dP 
dx 


h/ 


Eg = — 


/ {gop - gon) -{y-i) 7 ) dy 

1 ^ / 


i / {gop + gon)dy+ / (gop _go„) dy 

0 7 I 7 I 


(F3) 


where p, is the dynamic viscosity of the liquid wa¬ 
ter and A is the electrical conductivity of MgCl2 so¬ 
lution obtained from [s^. The sreaming currents per 
unit width from the MD simulations is calculated using 


Ii = Pel (y) ui (y) dy, where the variables pei (y) and 
Ml (y) represent local charge density and velocity calcu¬ 
lated from the MD simulations through post processing. 
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